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ABSTRACT 

The Mihalas-Hummer-Dappen (MHD) equation of state is a part of the 
Opacity Project (OP), where it mainly provides ionization equihbria and level 
populations of a large number of astrophysically relevant species. Its basic 
concept is the idea of perturbed atomic and ionic states. At high densities, 
when many-body effects become dominant, the concept of perturbed atoms 
loses its sense. For that reason, the MHD equation of state was originally 
restricted to the plasma of stellar envelopes, that is, to relatively moderate 
densities, which should not exceed p < 10^^ g cm~'^. However, helioseismological 
analysis has demonstrated that this restriction is much too conservative. 
The principal feature of the original [Hummer fc Mihalas (1988) paper is an 



expression for the destruction probability of a bound state (ground state or 
excited) of a species (atomic or ionic), linked to the mean electric microfield of 
the plasma. [Hummer fc Mihalas] (1988) assumed, for convenience, a simplified 



form of the Holtsmark microfield for randomly distributed ions. An improved 
MHD equation of state (Q-MHD) is introduced. It is based on a more realistic 
microfield distribution ( [Hooper[ 1966, 1968) that includes plasma correlations. 
Comparison with an alternative post-Holtsmark formalism (APEX) is made, 
and good agreement is shown. There is a clear signature of the choice of the 
microfield distribution in the adiabatic index 71, which makes it accessible 
to present-day helioseismological analysis. But since these thermodynamic 
effects of the microfield distribution are quite small, it also follows that the 
approximations chosen in the original MHD equation of state were reasonable. 
A particular feature of the original MHD papers was an explicit list of the 
adopted free energy and its first- and second-order analytical derivatives. The 
corresponding Q-MHD quantities are given in the Appendix. 
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1. Introduction 



The so-called MHD equation of state ( |Hummer fc Milialas] 1988; Milialas et al. 



|1988| ; Pappen et al. 1988| ; Pappen et al. 1987|) was developed as part of the international 
"Opacity Project" (OP; see [Seaton 1987, 1992, 1995| ; Perrington 1997] ). Its main purpose 



was to calculate the ionization degrees of all astrophysically relevant chemical elements in 
order to provide a crucial ingredient of the calculation of the radiative opacity of stellar 
interiors. The basic concept of the MHD equation of state was built on the idea of perturbed 
atomic and ionic states. At high densities, when many-body effects become dominant, the 
concept of perturbed atoms loses its sense. For that reason, the MHD equation of state 
was originally restricted to the plasma of stellar envelopes, that is, to relatively moderate 
densities, which should not exceed p < 10~^ g cm~'^. 

However, the MHD calculation of ionization equilibria was not only the necessary part 
of an opacity calculation. The same analytical and computational effort also allowed the 
computation of thermodynamic quantities to a high degree of accuracy and reliability. It 
turned out that for the purpose of thermodynamic calculations, the aforementioned density 
domain was much too conservative. For instance, [Christensen-Dalsgaard et al.| (1988) 
applied the MHD equation of state to models of the entire Sun in order to predict solar 
oscillation frequencies. MHD remained a reliable tool down to the solar center, where 
density is about 150 g cm~'^. So in spite of the original design of the MHD equation of state, 
the associated thermodynamic expressions have a broader domain of applicability, extending 
in particular to stellar cores. The reason is that in the deeper interior, the plasma becomes 
virtually fully ionized. Therefore, in practice, it does not matter that the condition for the 
legitimacy of the perturbation mechanism for bound species ( Pummer fc Mihala^ 1988) is 
not fulfilled, essentially because there are no bound species of the chemical elements that 
can be relevant for the equation of state (however, for opacity, bound states of less abundant 
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elements are relevant). Without bound species, the MHD equation of state falls back to 
an ideal-gas equation, enriched with Coulomb pressure and electron degeneracy. Coulomb 



pressure is only included to lowest order (the Debye-Hiickel approximation, see, e.g., pbeling 



et al. 197q ), yet it is now known that, somewhat fortuitously, the lack of higher-order 



Coulomb term has no significant consequence (for more details see [Christensen-Dalsgaard et 



al. 199(j| ). Specifically, it was shown by [Chabrier fc Baratfel (1997) that the MHD equation 



of state can be used for low-mass stars, at least down to 0.4 M© (see also [I'rampedach fc 
Uappen 1999|). 



Helioseismology has so far been putting the toughest observational constraints on the 
equation of state (for a recent review, see [Christensen-Dalsgaard et al. 199T; ). Given the 



importance of helioseismology as an experiment, we specifically review some of the most 
recent results in Section 3. It has emerged that the MHD equation of state clearly fares 
well, especially in comparisons with simpler formalisms, for instance the popular [Eggletoir 



Paulkner fc Flannery| (1973) equation of state (hereinafter EFF). However, there are still 
discrepancies between MHD and the inference from observations ([Christensen-Dalsgaard et 



|al. 1996| ). Some of these discrepancies were successfully removed by the OPAL equation 
of state (which is itself part of a major opacity project; [Rogers 198(j| ; [Iglesias fc Rogers 



1995; [Rogers et al. 19^ and references therein). However, recent results by [Basu et "al 



(1999) have indicated that in the outermost layers of the Sun {r / Rq > 0.97), the MHD 
equation of state appears to be a better match to the data than OPAL (see Section 3). 

The continual involvement of the MHD equation of state in current developments has 
led us to revisit the principal issues. The original [Hummer fc Mihalas[ (1988) paper derived 



an expression for the destruction probability of a bound state (ground state or excited) of a 
species (atomic or ionic). The probability was expressed as a function of the mean electric 
microfield. For convenience. Hummer fc Mihalas| (1988) assumed a Holtsmark microfield for 
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randomly distributed ions. Furthermore, to reduce the computational effort, they adopted 
a simplified approximate form of the Holtsmark function. 

It turned out that these simplifications have been the source of inadequacies of 
the MHD equation of state. In a real plasma, correlations occur, that is, the Coulomb 
interaction modifies the ion distribution. The result is that the microfield distribution peaks 
at lower values of the field strength than given by the Holtsmark distribution. Therefore, 
the probability with which a state for a specific atom or ion ceases to exist is reduced 
relative to the Holtsmark result, and thus finally, the occupation probabilities is higher 
than in the Holtsmark approximation. These discrepancies were demonstrated by [Iglesias 



fc Kogers| (1995). Therefore, the poor quality of the Holtsmark approximation is one of the 



reasons for the discrepancy between OP and OPAL opacities at higher densities. 

One might think that these differences would have only a small bearing on the equation 
of state ( [Rogers fc Iglesias| 1998). However, recent studies about the influence of details 



in the hydrogen partition function on thermodynamic quantities ( [Nayfonov fc Dappen 



1998; Pasu et aT] 1999) have revealed that the difference between various approximations 
of the microfield distribution is well within range of observational helioseismology (see 
Section 3). In view of these more stringent demands on the equation of state, we have 
improved the original MHD equation of state by including a post-Holtsmark microfield 
distribution to account for particle correlation. This microfield distribution function was 
derived by [Hoopei) (1966, 1968). Here, we present the resulting improved MHD equation 
of state. We name it "Q-MHD" after [Hummei]'s (1986) nomenclature (in which the 



microfield distribution was called P and its integral Q). Comparison with an alternative 
post-Holtsmark formalism (APEX; see [Iglesias et al. 1985| and references therein) is made. 



and good agreement is shown. Since one of the main features of the original MHD papers 
was an explicit list of the free energy and its first- and second-order analytical derivatives. 
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we list all the corresponding Q-MHD quantities in the Appendix. 



2. The MHD Equation of State 

The majority of the realistic equations of state that have appeared in the last 30 years 
are based on the so-called "free-energy minimization method". Before that, a common 
method was to use one expression {e.g., a modified Saha equation) for the ionization 
equilibria and another (not necessarily consistent) expression for the thermodynamic 
quantities. Such a procedure can lead to violations of Maxwell relations. The prevention 
is of course the use of a single thermodynamic potential, e.g., the free energy. Then, 
thermodynamic consistency is built in. However, since such a method requires considerable 
computing power, especially for multi-component plasmas, it could only become feasible 
since about 1960 (Parris 1959| ; [Harris et al. 196d| ). 



The free-energy minimization method uses statistical-mechanical models (for example, 
a partially degenerate electron gas, Debye-Hiickel theory for ionic species, hard-sphere 
atoms to simulate pressure ionization via configurational terms, quantum mechanical 
models of atoms in perturbed fields, etc.). It is a modular approach, that is, these models 
become the building blocks of a macroscopic free energy, which is expressed as a function 
of temperature T, volume V, and the particle numbers A^i, . . . , of the m components 
of the plasma. This model free energy is then minimized, subject to the stoichiometric 
constraints. The solution of this minimization problem then gives both the equilibrium 
concentrations and, if inserted in the free energy and its derivatives, the equation of state 
and the thermodynamic quantities. Obviously, this procedure automatically guarantees 
thermodynamic consistency. For obvious reason, this approach is called the "chemical 
picture" . Perturbed atoms must be introduced on a more or less ad-hoc basis to avoid the 
familiar divergence of internal partition functions (see, e.g., [Ebeling et al. 19T6D . 
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2.1. Internal Partition Functions 

The simplest way to modify atomic or ionic bound states to account for effects of the 
surroundings is by truncating the internal partition function at some maximum level that 
can be a function of temperature and density. Such a modification of bound states runs 
into the well-known technical problem that whenever the density passes through a critical 
value, for which a given bound state disappears into the continuum, the partition function 
changes discontinuously by the amount of the statistical weight Qi of the state. This is 
clearly unphysical and would lead to discontinuities and singularities in the free energy and 
its derivatives. 

One way of avoiding the problem of discontinuous jumps in the partition function is 
to assign "weights" or "occupation probabilities" to all bound states of all species. The 
internal partition function then becomes 

^jk = Yl^ijk 9ijkex.p [-Eijk/{kBT)] . (1) 

i 

Here, /cb is the Boltzmann constant, T temperature, Wijk the probability that the state i of 
ion j of species k still exists despite the plasma environment, and Qijk e'^^^''/^^^'^'' indicates 
the probability that this state is actually occupied in the system. Such an occupation 
probability formalism has several advantages: 

• The Wijk decrease continuously and monotonically as the strength of the relevant 
interaction increases. 

• States now fade out continuously with decreasing Wijk and thus assure continuity not 
only of the internal partition function but also of all material properties (pressure, internal 
energy, etc.) 

• The probabilistic interpretation of Wijk allows us to combine occupation probabilities 
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from statistically independent interactions. It is thus straightforward to allow for the 
simultaneous action of different mechanisms, accounting for several different species of 
perturbers by any one mechanism. Hence the method provides a scheme for treating 
partially ionized plasmas, and the limits of completely neutral or completely ionized gas are 
smoothly attained. 

• The Wijk can be made analytically different iable. In this way, MHD realized a reliable 
second-order numerical scheme in the free-energy minimization. 

• Finally, the Wijk formalism can be related naturally to line-broadening theory, which is 
important both for the interpretation of laboratory spectra and for opacity calculations. 

2.2. Statistical-Mechanical Consistency 

Although the Wijk can only be calculated a priori from some complementary interaction 
model, it must be stressed that one cannot introduce occupational probabilities into the 
internal partition function completely arbitrarily. 

Fermi (1924) showed that one can derive the Wijk from a free energy that depends 
explicitly on the individual occupation numbers Nijk] conversely it follows that the use of 
some heuristic w-ijk in the internal partition function implies the existence of an equivalent 
nonideal term in the free energy F (see [Fowler 1936|) . 

To illustrate that particular care must be taken when using Wijk-, consider a single- 
species perfect gas (that is, no dissociation or ionization), with a total of n particles, of 
which Hi are in state i (thus n = J^'^i = const). 

If / is defined by 
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Wi = e 



(df/dni)/(ksT) 



(2) 



then, to get a consistent free energy, the internal partition function 



■int 



[Eq. ©] 



must include the weights Wi of Eq. (Q). We denote by Z* such a particular internal 
partition function. For statistical-mechanical consistency, the free energy must contain a 
corresponding external part, which is the last term in the following equation for the free 



where V is the volume occupied by the gas, and A is the de Broglie wavelength of the 
particles of the species considered here. 

The last term in this equation has in fact been often ignored without justification. 
However, If the interaction term happens to be linear in the rij, then the last bracket 



the appropriately modified internal partition function Z* in the standard expression for 
F . Obviously, the bracket will not, in general, vanish for interaction terms that contain 
nonlinear combinations of the n^. However, in such cases it might be possible to find an 
astute simplified linear version of the interaction model which would reduce the problem to 
the previous case. 



energy of the gas ( [t'bwler 1936 ) 




(3) 




vanishes, and to get a consistent formulation it is sufficient merely to use 



2.3. MHD Occupational Probabilities 



2.3.1. Perturbations by Neutral Species 



For neutral perturbing species. Hummer fc IVlihala^ (1988) started out from a widely 
studied hard-sphere model ( [Fowler 1936| ) with each state in principle having its own 
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diameter. However, the simple binary interaction model is computationally prohibitive 
because it accounts for perturbations from all ions of all chemical species in all possible 
excited states. That implies thousands of the individual occupation numbers Nijk as 
independent variables in the free energy minimization. In addition, in this case the function 
/ is nonlinear. 

As an obvious first approximation, MHD considered the low- excitation limit (see 
Pummer fc Mihalas| 1988) in which it is assumed that essentially all perturbers encountered 



by an atom in an exited state reside in the ground state. 



(^Ufc)ne„tra1 = ^Xp 



j',k' 



(4) 



Here, the Vijk are the radii of state i of ion j of species k. Thus, the problem is reduced to 
total occupation numbers Nj^ = J2i ^ijk of all ionic species. And by eliminating the explicit 
appearance of Nijk the interaction / becomes linear in individual occupation numbers. 



/neutral = (|^) Nijk Y ^j'k'{rijk + Tijiyf . (5) 

^-^^ ^ i,j,k j',k' 

Therefore the term of the form 



i 

appearing in the multispecies generalization of Eq. (j^) vanishes identically, and the 
interaction only appears in the factors Wijk, that is, in Z*. 
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2.3.2. Perturbations by Charged Species 

In the case of charged particles, instead of trying to figure out / itself, MHD defines 
Wi directly, arguing that the presence of a plasma microfield destroys high-lying states by 
means of a series of Stark level mixing with higher lying states leading to the continuum. 

The basic idea is that for each bound state of every unperturbed ion (labeled by the 
indices i,j, k), there is a critical value of the electric field Fijk such that the state in question 
cannot exist if the field exceeds the critical value. Then probability that a given state does 
exist is simply the probabihty that the field strength is less than Fij^, i.e.. 



{Wijk 



yl^J charged 



ijk 



P{F)dF 



(6) 



where P{F) is the microfield distribution function (see Section 4). The choice of an 
appropriate plasma microfield P{F) is not straightforward. [Hummer fc Mihalas] (1988) 
have made the following choice (see section 4.2) of the resulting wijk, based on numerical 
comparisons with existing atomic physics calculations 



iJI^J charged 



exps-(^)16 



jk 



'^ijk Xij 



ijk 



j',k' 



(7) 



where Zjk denotes the charge of ion j of chemical species k (thus, zero for neutral particles) 
and the sum runs as before over all levels i of ions j of species k. Kijk is a quantum 
correction factor of those levels (see Eq. 4.70 of [Hummer fc Mihala^ 1988). Note that the 
interaction due to ^charged is automatically linear in the abundances Nij^ because this model 
of interaction does not depend on the internal excitation states of the perturbers. Therefore, 
there is again no term of the form of the last term in Eq. (^, and the interaction only 
appears in the factors Wijk- We stress that this property was of fundamental importance for 
the numerical feasibility of the MHD equation of state ([Hummer fc Mihalas| 1988). 
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Assuming statistically independent actions from neutral and charged perturbers, the 
joint occupation probability is the product of ("W^iifc) neutral (^^i^) charged- Formally, the 
neutral term could also be retained when extended charged particles interact with neutrals. 
However, because a hard-sphere description is highly implausible when one of the particles 
is charged, the MHD equation of state restricts the use of (w^jjA;)neutrai mutually neutral 
species (jk) only. 



3. Helioseismic Equation-of-State Diagnosis 

In the solar convection zone, helioseismology presents an opportunity to isolate 
the question of the equation of state from opacity and nuclear reaction rates, since the 
stratification is essentially adiabatic and thus determined by thermodynamics (see, e.g., 
Christensen-Dalsgaard fc Dappen 1992| ). Accurate analysis of the observations requires use 



of the full, nonasymptotic behavior of the oscillations. Figure shows a typical result of a 
numerical inversion ( [Basu fc (Jhristensen-Dalsgaard 1997| ). It shows the relative difference 
(in the sense Sun - model) between the squared sound speed obtained from inversion of 
oscillation data and that of a two standard models. A perfect model would lie on the zero 
line. The two models are in all respects identical except that they use a different equation 
of state (MHD and OPAL, respectively, cf. Section 1). Simplifying for the present purpose, 
we can look at results such as shown in Figure |l| as the data of helioseismology, disregarding 
how they were obtained from solar oscillation frequencies. 

EDITOR: PLACE FIGURE |l| HERE. 



The most important result of the earlier helioseismic equation-of-state analyses 
( Uhristensen-Dalsgaard et al. 199BD was that it is essential to include the leading Coulomb 
correction (the Debye-Hiickel term) to ideal-gas thermodynamics. Under solar conditions. 
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the size of the relative Coulomb pressure correction is largest in the outer part of the 
convection zone (about -8 percent) and it has another local maximum in the core (about -1 
percent). Figure |l], however, hides the fact that the Coulomb correction is the most 
important one, since both MHD and OPAL already contain it. Note that in Figure |l| 
the most significant information about the equation of state regards the convection 
zone (r > 0.71Rq), since beneath it, one cannot disentangle the influence from the 
equation of state from other effects. Figure Q contains the evidence that in the region 
O.9Oi?0 < r < O.97i?0, OPAL is a better fit to reality than MHD. 

Two very recent inversions have had further implications for the equation of state. 
First, the strong constraints from helioseismology now force us to include relativistic effects 
of electrons ( [Elliott fc Kosovichev 1998]) . Neither MHD nor OPAL have so far included 
relativistic effects, unlike the earlier and simpler EFF (Section 1) and its more recent sibling 
SIREFF ( Puzik fc Swenson 1997|) . Second, there are indications that for r > O.97i?0, 



MHD is favored over OPAL. This is the result of an apparent helioseismic confirmation 
of the thermodynamic effect of the excited states in hydrogen, treated according to MHD 
[Eq. (4), (6)]. The effect itself was demonstrated in a theoretical study by [Nayfonov 



|fc Dappen] (1998), which revealed interesting features beyond the Coulomb correction 
approximation. They are related to excited states and their treatment in the equation 
of state. The MHD equation of state with its specific, density-dependent occupation 
probabilities (Section 2.3) is causing a characteristic "wiggle" in the thermodynamic 
quantities, most prominently in xp = {dlnp/dln p)t, but equally present in the other 
thermodynamic quantities. 



EDITOR: PLACE FIGURE | HERE. 
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For convenience, in Figure 0, we recall the result for the helioseismically relevant 71 of 
a hydrogen-only plasma. Temperatures and densities are taken from a solar model. The 
appropriate density is implied but not shown (for more details and different thermodynamic 
variables and chemical compositions, see [Nayfonov fc Dappen| 1998 and Sections 6.2, 6.3). 



Five cases were considered: (i) MHD [standard MHD occupation probabilities: Eq. (4), (6)], 
(ii) MHDgs [standard MHD internal partition function of hydrogen but truncated to the 
ground state (GS) term], (iii) OPAL: OPAL tables [version of November 1996 of [Rogers et| 
|al. 1996| (1996)], (iv) MHDpL [MHD internal partition function of hydrogen, but replaced 
by the Planck-Larkin partition function Zpi = J2T Onl^^Pi—En/k^T) + E^/k^T — 1] 
( [Rogers 1986| )], (v) MHDpl,gs [MHDpL truncated to the ground state term]. The effect 



of the inclusion of the excited states in the internal partition function is manifest in the 
differences between MHD and MHDgs, and between MHDpL and MHDpl,gs, respectively. 
The effect of different occupation probabilities of ground and excited states shows up in 
the difference between MHD and MHDpL, and between MHDgs and MHDpl,gs- It was 
found that the presence of excited states is crucial. Also, the wiggle was demonstrated to 
be a genuine neutral-hydrogen effect despite the fact that most hydrogen is already ionized. 
This qualitative picture does not change when helium is added [Nayfonov fc Dappenj (1998). 

It seems that this effect of excited states can be observed in the Sun. Figure ^ shows 
the result of the inversion by jBasu et al] (1999), based on the solar oscillation frequencies 
obtained from the SOI/MDI instrument on board the SOHO spacecraft during its first 144 
days in operation ( [Rhodes et al.[ 1997). Previous sound-speed and 71 inversions {e.g., [Elliott 



[fc Kosovichev 19^) were still indirect, that is, they were giving the difference between solar 
models and the Sun without separating the change in structure due to the equation-of-state 
contribution. The inversion by [Basu et "aLl (1999) was done for the so-called "intrinsic" 



7i difference between solar models and the Sun. For convenience, we repeat the principal 
result in Figure which shows the difference between the Sun and various calibrated 
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solar models. The models alternate between two equations of state (MHD, OPAL), three 
different values for the solar radius, and two formalisms for convection [MLT: standard 
mixing length theory, CM: [Canuto fc Mazitelli] (1991) formalism]. The models are specified 
in Table 1, where in addition the calibrated values of the surface helium abundance Yg and 
the depth of the convection zone r^z are listed. All models assume gravitational settling of 
helium and heavy elements, an effect that is now part of the so-called standard solar model 
( lUhristensen-Dalsgaard et al. 1996| ). 



EDITOR: PLACE FIGURE | HERE. 



In contrast to earlier results for the intrinsic 71 difference ( [Basu fc Christensen- 
Dalsgaard 1997| ), which had uniform resolution throughout the Sun, the new study focuses 
on the 20% uppermost layers. It appears that in the top layers, the MHD models give a 
more accurate description of the Sun than the OPAL models. Since the difference in 71 
between MHD and OPAL is the wiggle of Figure ^ the observed preference of the MHD 
model in the upper region could indicate a validation of an MHD-like treatment of the 
exited states. Fig. |^ confirms the aforementioned earlier findings that below the wiggle 
region, OPAL fares better than MHD. In that region, the Planck-Larkin partition function 
of OPAL appears to be the better choice. 

The results in favor of MHD in the upper part of the Sun (r > 0.97Rq) could be 
due to the different implementations of many-body interactions in the two formalisms. 
Since density decreases in the upper part, OPAL by its nature of a systematic expansion, 
inevitably becomes itself more accurate; but MHD might, by its heuristic approach, have 
incorporated even finer, higher-order effects. 

Let us add a word of caution, though. It could appear tempting to produce a 
"combined" solar equation of state, with MHD for the top part and OPAL for the lower 
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part. However, such a hybrid solution is fraught with danger. For instance, it is known that 
patching together equations of state can introduce spurious effects ( [Dappen et al. 1993|) . 
It seems that the right way is to improve MHD and OPAL in parallel and independently, 
guided by the progress of helioseismology. The aim of the present study is to improve MHD. 



4. Microfield Distributions 

One of the mechanisms that affects the bound states of a radiator (atom or complex 
ion immersed in a plasma) is the electric field arising from the charged particles in its 
environment. 

To be able to calculate the internal partition function it, therefore, is necessary to 
estimate Enim and to calculate a microfield distribution. The actual electric field at a 
radiator in plasma is fluctuating in time. From the time-scale point of view, it is convenient 
to split the electric field into two components: a high-frequency part (with respect to the 
time scale of particle collisions) and a low-frequency component which varies on a time 
scale much longer than the orbital period of the bound state considered. Because of their 
high velocities, free electrons are considered as perturbers in the high-frequency part of the 
microfield distribution only. The distribution of the low-frequency component is calculated 
by considering a gas of ion perturbers. As was shown in ( [Hummer fc Mihalas| 1988), the 



low-frequency component dominates microfield effects at least for the conditions of stellar 
envelopes. In the following we discuss only this component. 



4.1. Holtsmark Microfield Distribution 



Unsold (1948) suggested a model for a microfield distribution that allowed for a 
simple analytical formulation. It contained the nearest neighbor (NN) approximation and. 
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therefore, was better suited to study the hmit of strong fields. The long range of Coulomb 
interactions in plasmas, however, renders the nearest neighbor approximation inadequate. 
The next logical step in the development of the microfield distributions should involve the 
interaction of the radiator with all ions in plasma. For a case of pure hydrogen plasma this 
problem was worked out by |Holtsmark| (1919). Two serious shortcomings of the Holtsmark 
distribution are: 

• its limitation to neutral radiators and 

• the absence of correlations between the charged plasma perturbers. 

The original MHD formalism ( [Hummer fc Mihala^ 1988) modified the Holtsmark 



distribution by making plausible but non-rigorous correction for the effects of ions other 
than protons. It considers hydrogenic radiators in a plasma for which the microfield 
perturbations by a variety of ionic species of charge Zp, p = 1,2, are dominant. The 
resultant occupation probability of an electron in level i = (n, /, m) of a hydrogenic potential 
of charge Za is given by (see Eq. 4.68a of [Hummer fc Mihalas| 1988) 



where 



«Aj / ^''"'e \ / "'ion \ 



^ZaalK 3 



Hp 



Q{x)= r PhWI^ . (9) 

JO 

Ki is the quantal Stark-ionization correction factor, Xi is the ionization potential of the 
level z, is the Bohr radius, Uf. is the electron density, rtion is the density of all ions, and 
Ph{(3) is the Holtsmark distribution function 
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Ph{P) = (v) l^dyexp{-y'/')ysm(3y. (10) 

4.2. MHD Microfield Distribution 

In the form of Eq. (|^) the calculations of plasmas with dozens of ions and hundreds of 
bound states are too cumbersome, especially if one wants to code all the derivatives in the 
analytical form. Instead, MHD relied on numerical experiments in search for analytical fits 
that would mimic the Holtsmark distribution function. It turned out that a good starting 
point was the |Unsold| (1948) occupation probability w'^, given by 



It might be worth noting here that because Xi V''^^^ even putting the arbitrary factor 
of 2 into exponent of Eq. (11) changes the cutoff quantum number by only a factor of 
2^/^ = 1.12, which is typically quite negligible. 

This simple expression is indeed a good fit to the Holtsmark distribution, if an 
additional more-or-less ad hoc factor 2 is put in the argument of Eq.(^. Tests showed 
( [Hummer fc Mihalas] 1988) that for strongly bound states with wf > 0.1 this trick actually 



leads to a good fit. However, for wf < 0.1 this fit decreases much more sharply than does 
Eq.(|D; for many applications this is not a matter of serious concern inasmuch as the basic 
physical effect - that such levels are for all practical purposes destroyed - is achieved. 

The adopted form of in MHD equation of state reflects the interactions of the 
radiator with all ions in the plasma, but the approximation does not take into account 
the correlations of the charged perturbers. This effect is already a problem because of 
the choice of the Holtsmark distribution function, not only one of the replacement of the 
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Holtsmark distribution by a modified Unsold expression (|lglesias fc Rogers] 1995). 



4.3. Microfield Distributions with Correlation Effects 

4.3.1. APEX Microfield Distribution 

One of the alternative methods to calculate the microfield distribution for plasmas of 
general nature, which include particle correlations, is the APEX ("adjustable-parameter 
exponential" ) approximation ( see [iglesias et al. 19"55| and references therein) based on a 



so-called "independent quasiparticle model". 

Originally developed as a phenomenological method, APEX was shown to be in 
agreement with a rigorous theoretical procedure (Putty et al. 198"5| ) related to a Baranger- 
Mozer series type of analysis. APEX has been shown to agree quite well with computer 
Monte-Carlo simulations for both high- and low-frequency component distributions. 
Especially well-suited for high-Z plasmas ( [Iglesias et al. 198^) , it treats the low-frequency 
component by considering a gas of ions interacting through electron screened potentials, 
which is a way to include static contributions from both ions and electrons into account. 
The success of APEX can be attributed to a fact that for small fields the contributions 
from many ions are important and these are well characterized by the second moment of 
the microfield distribution that is exactly included in APEX. 

The correlation effects are incorporated through the radial distribution functions of 
the ionic perturbers calculated in the hyper-netted-chain approximation generalized to 
multi-component plasmas ( [Rogers 1980| ). 



- 21 - 



4-3.2. Q-fit Micro field Distribution 

An alternative improvement over the Holtsmark distribution function, which accounted 
for particle correlations yet remained sufficiently simple so that it could be expressed by 
analytical fits suited to an equation-of-state program, was introduced by one of the authors 
(D. Hummer) in 1990. The numerical fits were based on work by Pooper| (1966, 1968). 
Although until the present work these fits have never been included in equation of state 
calculations, they were successfully used in an non-LTE investigation of the stellar fiux near 
the series limits in hot stars ( jHubeny, Hummer fc Lanz 19941 ). 



The |Hooper| (1966, 1968) microfield distribution was derived under the following 



assumptions: 

• the perturbing ions and electrons are in equilibrium with the same kinetic temperature, 

• while multiply charged radiators were allowed, all of the perturbing ions were singly 
charged. 

While these assumptions are obviously too stringent for studies of general plasmas and 
even laser-generated plasmas (see P?ighe fc Hooper 1977|) , they should be quite reasonable 
for the plasma of the interior of more or less normal stars having typical astrophysical 
compositions. In this case, one can show that the two assumptions of the Hooper microfield 
are indeed satisfied. First, the interior of normal stars is in thermodynamic equilibrium. 
Electrons and ions have therefore the same temperature. Note that stellar atmospheres are 
of course not in thermal equilibrium, but any stellar equation of state of the type discussed 
in this article is inadequate for them, and special nonlocal formalisms are needed. Second, 
for not-too-evolved stars, the typical chemical composition of the universe prevails, which 
has about 90% hydrogen by number. For such plasmas, the vast majority of perturbing 
ions are protons, i.e., singly charged. 
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When plasma correlation effects are important, the microfield distribution function 
depends on two additional parameter, the radiator charge and the correlation parameter 



-J'^^^d) = , (12) 



where Up is in cm ^ and T is in K. 



The physical meaning of the parameter a is clear from Eq. (0), namely a = rj where 
is the number of ions in a Debye sphere of radius ro 



In the low-frequency approximation only electrons contribute to shielding [Eq. (pTSl)] 
and the overall charge neutrality of plasma was used to obtain the last form of Eq. ([12|). 

Another way to look at a is to realize that it can be expressed in terms of the electron 
coupling parameter Fei 

where Oe is the electron sphere radius defined by ^rigal = 1. 

The appropriate microfield distribution function W{l3] Zr, a) has been derived 
by |Hooper| (1966, 1968). Using a substantially modified version of Hooper's code, Thomas 
Schoning of the University of Munich Observatory and Hummer computed two-dimensional 
fits in (3 and a to W{f3; Zr, a) for = 0, . . . , 5 and a < 0.8. From these fits Hummer 
evaluated the function defined by Eq. (^ 



Q{P; Zr, a) = f W{(3'- Zr, a)d(3' 
Jo 



(15) 
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for the given ranges of the parameters and a. It was possible to get a reasonably good 
analytical fit to obtained data by using just two parameters: radiative charge Z^. and a 
correlation parameter a. 

The adopted form of the fit Q{/3,Zr,a) is 



where 



The coefficients C\ and C2 depend on a and Z^. through the forms: 

Cx{Zr, a) = Pi [X + P^ZrC?\ , (18) 

and 

C2 = P2X, (19) 

where 

X ^ {1.0 + Psaf\ (20) 
and the optimum values of the parameters seem to be: 



Pi 
P2 
P3 
Pa 
P5 



0.1402 
0.1285 
1.0 
3.15 
4.0 . 



(21) 
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The chosen forms of the fits are constrained to lie between zero and 1 for all values of 
a and and to give the correct functional form for both small and large limiting values 
of (3. The presented fits have been obtained for data from a = (no correlation, which 
corresponds to the Holtsmark limit) to a = 0.8. For a=Zr=0, the fit is accurate to within 
±2%. Except for very small values of /5, the fit is accurate to within 10%, except for large 
a and Z^., where it reaches 26% in the worst case ( for very small (3). However, as Q{I3) is 
very small there ( on the order of 10~^ ), it shouldn't matter. Figures ^ and |^ demonstrate 
the magnitude of the correlation effects for a neutral as well as for a singly charged radiator. 

In the following we are going to ignore the temperature and density dependence of the 
coefficients Ci,C2 in the derivatives of the Q{(3). For solar conditions, the maximum relative 
error of this approximation is about 10^^-10^^ and, therefore, the approximation is quite 
sufficient for our fits, which in any case have a few percent errors themselves. Henceforth, 



df _ 2C,P\l + \C^^^I') 
d(3 (1 + C2P'/^f 



(22) 



4 {l + C2f3^/^f 



(23) 




(24) 




(25) 



5. Q-MHD Equation of State 
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5.1. Q-fit Occupation Probability 

The post-Holtsmark microfield distribution given by tlie Q-fit can be used to upgrade 
tlie original MHD to the Q-MHD equation of state. Recall (Section 2.3) that the MHD-style 
occupation probability is given by a product of neutral and charged parts 

Wijk = w-jkWijk ■ (26) 

The neutral part wfjj., corresponding to pressure ionization mechanism, has the same 
form as in the original MHD representation (see pummer fc Mihalas] 1988). The charged 



part wfjfc, determined by the microfield distribution, is given by 



<k = QW^^k) , (27) 



_ / 3 K^jkXlk fn^onV/' 1 , 

Utt) 4{Z,, + l)eAnJ nT ' ^ 

where Xijk is the ionization potential of level i of ion j of chemical species fc, is the 
net charge on ion j of species k ( zero for neutral particle), rijo^ is the total density 
of ionic perturbers in a system and is again the electron density. is the quantal 
Stark-ionization correction factor of MHD theory (see [Hummer fc Mihalas| 1988) 



Kn = 1 , n <3, 

Kn = n>3. (29) 

3 {n + 1) 

In the following we adopt an important approximation {nion/n^)^^^ ^ 1. As been 
already discussed in [Hummer fc Mihalas] (1988) the error of this assumption for stellar 
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plasmas of normal chemical compositions never exceeds a few percent due to scarcity of 
high-Z species. This minor approximation simplifies the analytical derivatives (Appendix) 
significantly. 

Therefore, Equation (^8[) can be rewritten as 



5.2. Q-MHD Free Energy 

Once the new Q-fit occupation probabilities are known, they merely have to be put in 
place of the original MHD occupation probabilities. To allow detailed comparisons with the 
MHD expressions ( |Mihalas et al. 1988 ; Pappen et al. 1988 ), we list the Q-fit analogs in the 



Appendix. Calculation of the ionization equilibria and thermodynamic quantities is done 
as in MHD, and the result is the Q-MHD equation of state. 



6. Results and Discussion 

6.1. Comparisons of Q-fit and APEX Microfield Distributions 

In this section the microfield distributions in hydrogen plasmas are presented for 
different values of the plasma coupling parameter F (which is the ionic analog to Fe defined 
by Eq. (|1^. More precisely. 




where Oav is the average ion sphere radius defined by (47r/3) Z^ion^j'^av = 1- both Q-MHD 
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and APEX formalisms only singly-charged perturbers considered. To demonstrate the effect 
of a radiator charge, the case of a neutral {Z^ = 0) radiator (see Figs. | and ^ is presented 
alongside with the singly- charged {Zr = 1) case (Figs. calculated for a several values of 
coupling parameter F. 

EDITOR: PLACE FIGURE | HERE. 
EDITOR: PLACE FIGURE | HERE. 
EDITOR: PLACE FIGURE | HERE. 

EDITOR: PLACE FIGURE HERE. 

The resultant function Q [see Eq. (^] for these cases demonstrates a clear dependence 
on a radiator charge as well as on a value of F. 

EDITOR: PLACE FIGURE | HERE. 
EDITOR: PLACE FIGURE HERE. 



EDITOR: PLACE FIGURE 10 HERE 



EDITOR: PLACE FIGURE 



11 



HERE. 
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6.2. Comparison of the Equation of State for a Hydrogen-only Plasma 

This section deals with equation-of-state calculations for hydrogen-only plasmas. Four 
different equations of state have been calculated for a set of fictitious solar temperatures 
and densities, given in Fig. The densities were chosen by [Nayfonov fc Dappen] (1998) 



to simulate solar pressure at the given temperature for a hydrogen-only plasma. In the 
following we will denote by "H-only solar track" the set of temperatures and densities of 
Figure O. Figure O demonstrates a coupling parameter F as estimated for the conditions 



of the H-only solar track. 



EDITOR: PLACE FIGURE 12 HERE. 



EDITOR: PLACE FICURE 13 HERE. 



The four equations of state (EOS) are: (i) regular MHD EOS (solid lines in all of the 
following graphs); (ii) MHD EOS with a true Holtsmark microfield distribution function 
(dotted-dashed lines); (iii) Q-MHD EOS, as described in Section 5 of this paper (dashed 
lines) and, finally, (iv) OPAL EOS (dotted lines). The results expand the previous study 
( [Nayfonov fc Dappen] 1998) on the effects of internal partition functions, and in as much 
as the studies overlap, they agree with each other. Among the quantities considered, that 
is, Xp = {dlnp/dln p)j., xt = {d\np/d\nT)^, and 71 = {d\np/d\n p)^, only Xp reveals 
differences already in the absolute plot. However, the differences are present, and of the 
same order of magnitude, also in the other two quantities, in the same way as found by 
Nayfonov fc Dappen] (1998). Figures |l^ and |l^ show absolute values of Xp and relative 



difference (with respect to OPAL), respectively. All three MHD-type EOS {i.e., MHD, 
MHD with true Holtsmark, and Q-MHD) demonstrate the characteristic wiggle discovered 
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in the previous study (see also Fig. ^. The wiggle is located in the temperature zone 
4.5 < logT < 5.3. Given the hydrogen-only plasma, it is obviously a pure hydrogen 
phenomenon. However, as the location of the 10% to 90% ionization zone demonstrates 



(Fig. |T5|), the wiggle is caused by the last remaining neutral hydrogen atoms, fighting 



against pressure ionization in a region of near full ionization. 



EDITOR: PLACE FIGURE 14 HERE. 



EDITOR: PLACE FIGURE 15 HERE. 



The figures also reveal a close agreement between regular MHD and MHD with a true 
Holtsmark distribution function, which gives a post factum justification of the choice made 
by the authors of MHD in 1988. It is also clear that Q-MHD seems to be in a better overall 
agreement with OPAL than other two equations of state. 



6.3. Comparison of the Equation of State for a Hydrogen-Helium Mixture 

To study a more realistic picture of solar plasmas, calculations similar to ones described 
in the previous section were carried out for a hydrogen-helium mixture along a real 
solar profile of density and temperature (model S of [Christensen-Dalsgaard et al. 199^ . 



Hydrogen constitutes 74% of this mixture by mass, making it very similar to a solar 
composition of the standard solar model. While those calculations exclude heavier elements 
from the analysis, the obtained results actually give a strong indication that a treatment 
of heavier elements is necessary to improve our current equation of state models ( |Dappen 



let al. 199^ ). The plots of absolute values of Xp and of relative differences (with respect to 



OPAL) between the same four equations of state are given by Figures |T6| and [T^. 
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EDITOR: PLACE FIGURE |T6| HERE 



EDITOR: PLACE FIGURE ITTI HERE 



EDITOR: PLACE FIGURE |T8| HERE. 



Figure |T8| reveals that the signature of the hydrogenic internal partition function is 
present in other thermodynamic quantities, such as 71, as well. Again, one can argue that 
Q-MHD agrees with OPAL better than the other models and that the differences between 



the usual MHD and MHD with a true Holtsmark distribution remain small. Fig. \18 
shows that Q-MHD leaves MHD and approaches OPAL especially in the temperature 
range 5.1 < logT < 5.4. From the comparison with the helioseismic experiment (Fig. 0), 
one realizes that for at least the upper part of this temperature range, OPAL is a better 
equation of state than MHD. Therefore, Q-MHD corrects MHD in the right direction. 

7. Conclusion 



Upgrading the MHD equation of state to include realistic microfield distributions 
beyond the Holtsmark approximation has confirmed the significant changes in the 
occupation numbers for atomic and ionic states, as was expected by |lglesias fc Rogers 



(1995). These changes in the occupation numbers and the associated shifts in the ionization 
balances are widely assumed responsible, among other, for the discrepancies between OP 
and OPAL opacities under the temperatures and densities of the solar center (see Gong 



et al. 19981) . For condit 10ns of envelopes of more massive stars, however, the two opacity 



calculations agree very well. Our study is a systematic comparison of the impact on 
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occupation probabilities of the original MHD microfield, the proper Holtsmark microfield, 
and the APEX distribution. 



As far as thermodynamic properties are concerned, Iglesias fc Rogers] (1995) and also 



Rogers fc Iglesias| (1998) believed that effects of different microfield distributions would 



be rather negligible. But encouraged by recent progress on the influence of excited states 
in hydrogen on thermodynamic quantities (Section 3), we have found a clear signature 
of the microfield distribution, easily within reach of helioseismological accuracy (compare 



Figures ^ and |18D . This confirms once more that solar observations constrain formalisms 



used to describe the physics of atoms and compound ions immersed in a plasma. 

Since these thermodynamic effects are nonetheless quite small and are, even for 
helioseismological accuracy, only relevant at some selected locations of the Sun, our results 
also show that for most applications of stellar structure, the approximations chosen in the 
original MHD equation of state are reasonable. In the particular cases where they are not 
reasonable, such as in helioseismological studies of the zones of partial ionization in the 
Sun, the more accurate Q-MHD equation of state, which is based on a realistic microfield 
distribution, is more accurate. 
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A. Appendix: Analytical Expressions for the Free Energy and Its Derivatives 

with Q-fit Occupation Probabilities 



We follow the MHD notation (Mihalas et al. 1988; Dappen et al. 1988) very closely. 



The term affected by the new microfield is the free energy of bound systems of species s 



F2 = J2Ns {Els - ksT log Zs) , (Al) 
with the internal partition function evaluated for energies relative to the ground state, that 



Zs = J2 ^^s9is exp ^ . (A2) 

Here, Ai^j^ = Eis — Eig. In what follows, A,/i, u label neutral, q charged particles. We shall 
also treat (Np/Ne) as a constant (see the discussion leading to Eq. (0). Then, (1) F2 is 
linear in Nj, and (2), Zg does not depend on individual ions Nj. Then 



||^^.-.er,nZ„-.eTg||| (A3) 



^ = E,^-k^T\nZ, (A4) 
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dV 



-AEi,/kBT 



(A25) 



d'Z,, 1 



dTdV k^T^V 



(A26) 
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3 
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47r 



-AEi,/kBT 
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Table 1: Properties of the solar models of Fig. |3|. See tex t. 

Model EOS Radius Convective Yg Vcz/ Rq 

Mm Flux 

Ml MHD 695.78 CM 0.2472 0.7145 

M2 MHD 695.99 CM 0.2472 0.7146 

M3 MHD 695.51 CM 0.2472 0.7145 

M4 MHD 695.78 MLT 0.2472 0.7146 

M5 OPAL 695.78 CM 0.2465 0.7134 

M6 OPAL 695.99 CM 0.2465 0.7135 

M7 OPAL 695.51 CM 0.2466 0.7133 

M8 OPAL 695.78 MLT 0.2465 0.7135 
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Fig. 1. — Difference between squared sound speed from inversion of oscillation data and that 
of a standard model based on the MHD (circles) and OPAL (triangles) equation of state 
(Figure provided by S. Basu). 
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Fig. 2. — Left Panel: absolute values of 71 for solar temperatures and densities of a hydrogen- 
only plasma. Linestyles: MHD - asterisks, MHDqs - dashed hues, MHDpL - dotted-dashed 
lines, MHDpL,GS ~ dotted lines, and OPAL - solid lines. See text for the definitions of the 
different MHD versions. Right Panel: relative differences with respect to MHDgs, in the 
sense (71 — 71 [MHDgs] )/7i[MHDgs]), using the same hne styles as in the left panel. The 
horizontal solid zero line, representing MHDgs, is also shown. 
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Fig. 3. — Relative difference between 71 obtained from an inversion of lielioseismological 
data and 71 of tlie solar models listed in Table I, in the sense "Sun - model". Only the 
intrinsic difference in 71 is shown, that is, the part of the difference due to the equation of 
state (Figure provided by S. Basu). 
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Fig. 4. — Microfield distribution from Q-MHD model in a case of neutral perturber in 
hydrogen plasma for different values of coupling parameter: F = 0.10 (solid line), F = 0.25 
(dashed line), F = 1.0 (dotted-dashed line). 
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Fig. 5. — Microfield distribution from APEX model in a case of neutral perturber in hydrogen 
plasma for different values of coupling parameter: F = 0.10 (solid line), F = 0.25 (dashed 
line), F = 1.0 (dotted-dashed line). 
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Fig. 6. — Same as Fig.l in a case oi Zr = I perturber. 
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Fig. 7. — Same as Fig. 2 in a case oi Zr = I perturber. 




-48- 



Fig. 9. — Function Q for a case presented in Fig. 2 
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Fig. 12. — "H-only solar track" : for given temperature, density of a hydrogen-only plasma 
is chosen such that pressure corresponds to solar pressure. 
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Fig. 13. — Estimated coupling parameter T along the H-only solar track 
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Fig. 14. — Xp foi" ^ case of hydrogen-only plasma calculated for 4 different models: standard 
MHD (solid line), MHD with Holtzmark (dotted-dashed line), Q-MHD (dashed hne) and 
OPAL (dotted line). For details see text. 
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Fig. 15. — Relative differences in xp with respect to OPAL for a case of liydrogen-only 



plasma. Linestyles are the same as in Fig. |Tj. The ionization zone of hydrogen (ionization 
degree between 10% and 90%) is indicated. For details see text. 
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Fig. 16. — The same as Fig. for a case of a hydrogen-helium mixture (74% to 26% by 
mass, respectively). For details see text. 
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Fig. 17. — The same as Fig. |T5| for the hydrogen-hehum mixture of Fig. Besides the 
hydrogen ionization zone, the two ionization zones of hehum (ionization degree between 10% 
and 90%) are also indicated. For details see text. 
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Fig. 18. — Relative differences in the adiabatic index 71 with respect to OPAL for the 



hydrogen-hehum mixture of Fig. |TJ. Ionization degrees are as in Fig. |T^, linestyles as in 
Fig. |l4l For details see text. 



